library(mvtnorm)
library(foreign)

###############
#  Procedures #
###############
 

# IV probability test based on Bennett (2009) and Chen and Schroeter (2012)
ivtestprob<-function(y,d,z,n_boot1=999, n_boot2=999, gridval=NULL, subsets=5){
n=length(y)
p1_1<-mean(d[z==1])
p1_0<-mean(d[z==0])
p0_1<-1-mean(d[z==1])
p0_0<-1-mean(d[z==0])
z1<-z[d==1]
z0<-z[d==0]
y1<-y[d==1]
y0<-y[d==0]
y11=y1[z1==1]
y01=y1[z1==0]
y10=y0[z0==1]
y00=y0[z0==0]

if ((is.null(gridval)*1)==1){
gridval=seq(from=min(y),to=max(y),length.out = subsets+1)
grid=subsets
}
if ((is.null(gridval)*1)==0){
grid=length(gridval)-1
}
p11<-c()
p01<-c()
p10<-c()
p00<-c()

gridval[length(gridval)]=gridval[length(gridval)]+0.00001
for (i in 1:(grid)){
p11<-c(p11,mean(y11>=gridval[i] & y11<gridval[i+1]))
p01<-c(p01,mean(y01>=gridval[i] & y01<gridval[i+1]))
p10<-c(p10,mean(y10>=gridval[i] & y10<gridval[i+1]))
p00<-c(p00,mean(y00>=gridval[i] & y00<gridval[i+1]))
}
lengthp=length(p11)

q<-p1_0/p1_1
r<-p0_1/p0_0

p11min<-(p11-(1-q))/q
p11max<-p11/q
p00min<-(p00-(1-r))/r
p00max<-p00/r


#first bootstrap 
Tn1b = c()
Tn2b = c()
Tn3b = c()
Tn4b = c()

while(length(Tn1b)<n_boot1*lengthp){
sboot<-sample(1:length(y),length(y),TRUE)
zb<-z[sboot]
db<-d[sboot]
yb<-y[sboot]

obsb=length(yb)
obs1b<-sum(db)
obs0b<-sum(1-db)

p1_1b<-mean(db[zb==1])
p1_0b<-mean(db[zb==0])
p0_1b<-1-mean(db[zb==1])
p0_0b<-1-mean(db[zb==0])

if ((p1_0b/p1_1b<=1) & (p0_1b/p0_0b<=1)){
z1b<-zb[db==1]
z0b<-zb[db==0]
y1b<-yb[db==1]
y0b<-yb[db==0]
y11b=y1b[z1b==1]
y01b=y1b[z1b==0]
y10b=y0b[z0b==1]
y00b=y0b[z0b==0]

p11b<-c()
p01b<-c()
p10b<-c()
p00b<-c()
for (i in 1:(grid)){
p11b<-c(p11b,mean(y11b>=gridval[i] & y11b<gridval[i+1]))
p01b<-c(p01b,mean(y01b>=gridval[i] & y01b<gridval[i+1]))
p10b<-c(p10b,mean(y10b>=gridval[i] & y10b<gridval[i+1]))
p00b<-c(p00b,mean(y00b>=gridval[i] & y00b<gridval[i+1]))
}

qb<-p1_0b/p1_1b
rb<-p0_1b/p0_0b
p11minb<-(p11b-(1-qb))/qb
p11maxb<-p11b/qb
p00minb<-(p00b-(1-rb))/rb
p00maxb<-p00b/rb

Tn1b = cbind(Tn1b,(p11minb-p01b)) 
Tn2b = cbind(Tn2b,(p01b-p11maxb)) 
Tn3b = cbind(Tn3b,(p00minb-p10b)) 
Tn4b = cbind(Tn4b,(p10b-p00maxb)) 
}
}

# variances
varTn1b=diag(var(t(Tn1b)))
varTn2b=diag(var(t(Tn2b)))
varTn3b=diag(var(t(Tn3b)))
varTn4b=diag(var(t(Tn4b)))

# create objects for partial and full recentering
FCbsample1 <-matrix(0,nrow=lengthp,ncol=n_boot2)
FCbsample2 <-matrix(0,nrow=lengthp,ncol=n_boot2)
FCbsample3 <-matrix(0,nrow=lengthp,ncol=n_boot2)
FCbsample4 <-matrix(0,nrow=lengthp,ncol=n_boot2)

PCbsample1 <-matrix(0,nrow=lengthp,ncol=n_boot2)
PCbsample2 <-matrix(0,nrow=lengthp,ncol=n_boot2)
PCbsample3 <-matrix(0,nrow=lengthp,ncol=n_boot2)
PCbsample4 <-matrix(0,nrow=lengthp,ncol=n_boot2)

FCp_vals1 <- matrix(0, nrow = n_boot2, ncol = lengthp)
FCp_vals2 <- matrix(0, nrow = n_boot2, ncol = lengthp)
FCp_vals3 <- matrix(0, nrow = n_boot2, ncol = lengthp)
FCp_vals4 <- matrix(0, nrow = n_boot2, ncol = lengthp)

PCp_vals1 <- matrix(0, nrow = n_boot2, ncol = lengthp)
PCp_vals2 <- matrix(0, nrow = n_boot2, ncol = lengthp)
PCp_vals3 <- matrix(0, nrow = n_boot2, ncol = lengthp)
PCp_vals4 <- matrix(0, nrow = n_boot2, ncol = lengthp)

#Fix value of sequence delta_n
delta_n = sqrt(2*log(log(n))/n)

# compute test statistics
Tn1 = (p11min-p01) 
Tn2 = (p01-p11max)
Tn3 = (p00min-p10) 
Tn4 = (p10-p00max)  
Tn=c(Tn1, Tn2, Tn3, Tn4) 

# compute recentering terms
maxVec1=( Tn1 >= (-delta_n*sqrt(varTn1b))  )*Tn1    + ( Tn1 <  (-delta_n*sqrt(varTn1b))  )* (-delta_n*sqrt(varTn1b) )
maxVec2=( Tn2 >= (-delta_n*sqrt(varTn2b))  )*Tn2   +( Tn2 < (-delta_n*sqrt(varTn2b))  )* (-delta_n*sqrt(varTn2b) )
maxVec3=( Tn3 >= (-delta_n*sqrt(varTn3b))  )*Tn3    + ( Tn3 <  (-delta_n*sqrt(varTn3b))  )* (-delta_n*sqrt(varTn3b) )
maxVec4=( Tn4 >= (-delta_n*sqrt(varTn4b))  )*Tn4   +( Tn4 < (-delta_n*sqrt(varTn4b))  )* (-delta_n*sqrt(varTn4b) )

Tn1 = sqrt(n)*Tn1 
Tn2 = sqrt(n)*Tn2
Tn3 = sqrt(n)*Tn3 
Tn4 = sqrt(n)*Tn4

Tn1b=sqrt(n)*Tn1b
Tn2b=sqrt(n)*Tn2b
Tn3b=sqrt(n)*Tn3b
Tn4b=sqrt(n)*Tn4b

# Generate bootstrap samples (using first bootstrap)- full recentering
FCBootstrap1<-matrix(Tn1b-rep(Tn1,n_boot1),lengthp,n_boot1)
FCBootstrap2<-matrix(Tn2b-rep(Tn2,n_boot1),lengthp,n_boot1)
FCBootstrap3<-matrix(Tn3b-rep(Tn3,n_boot1),lengthp,n_boot1)
FCBootstrap4<-matrix(Tn4b-rep(Tn4,n_boot1),lengthp,n_boot1)

# Generate bootstrap samples (using first bootstrap)- partial recentering
PCBootstrap1<-matrix(Tn1b-rep((sqrt(n)* maxVec1),n_boot1),lengthp,n_boot1)
PCBootstrap2<-matrix(Tn2b-rep((sqrt(n)* maxVec2),n_boot1),lengthp,n_boot1)
PCBootstrap3<-matrix(Tn3b-rep((sqrt(n)* maxVec3),n_boot1),lengthp,n_boot1)
PCBootstrap4<-matrix(Tn4b-rep((sqrt(n)* maxVec4),n_boot1),lengthp,n_boot1)

# Generate empirical distribution functions - full recentering
FC_ECDF1 <- apply(FCBootstrap1, 1, ecdf)
FC_ECDF2 <- apply(FCBootstrap2, 1, ecdf)
FC_ECDF3 <- apply(FCBootstrap3, 1, ecdf)
FC_ECDF4 <- apply(FCBootstrap4, 1, ecdf)

# Generate empirical distribution functions - partial recentering
PC_ECDF1 <- apply(PCBootstrap1, 1, ecdf)
PC_ECDF2 <- apply(PCBootstrap2, 1, ecdf)
PC_ECDF3 <- apply(PCBootstrap3, 1, ecdf)
PC_ECDF4 <- apply(PCBootstrap4, 1, ecdf)

# second bootstrap - draw observations from first bootstrap
bsample2 <- sample(n_boot1, n_boot2, replace=TRUE)

# draw samples of recentered test statistics out of first bootstrap - full recentering
FCbsample1 = FCBootstrap1[1:lengthp,bsample2]
FCbsample2 = FCBootstrap2[1:lengthp,bsample2]
FCbsample3 = FCBootstrap3[1:lengthp,bsample2]
FCbsample4 = FCBootstrap4[1:lengthp,bsample2]

# draw samples of recentered test statistics out of first bootstrap - partial recentering
PCbsample1 = PCBootstrap1[1:lengthp,bsample2]
PCbsample2 = PCBootstrap2[1:lengthp,bsample2]
PCbsample3 = PCBootstrap3[1:lengthp,bsample2]
PCbsample4 = PCBootstrap4[1:lengthp,bsample2]

# empirical distributions of fully recentered p-values

p_val_1 =1
p_val_2 =1
p_val_3 =1
p_val_4 =1

for (k in 1:lengthp) {
	FCp_vals1[1:n_boot2,k] = 1-FC_ECDF1[[k]](FCbsample1[k,1:n_boot2])	
	FCp_vals2[1:n_boot2,k] = 1-FC_ECDF2[[k]](FCbsample2[k,1:n_boot2])
	FCp_vals3[1:n_boot2,k] = 1-FC_ECDF3[[k]](FCbsample3[k,1:n_boot2])	
	FCp_vals4[1:n_boot2,k] = 1-FC_ECDF4[[k]](FCbsample4[k,1:n_boot2])
	FCp_vals=cbind(FCp_vals1,FCp_vals2,FCp_vals3,FCp_vals4)

# empirical distributions of partially recentered p-values
	PCp_vals1[1:n_boot2,k] = 1-FC_ECDF1[[k]](PCbsample1[k,1:n_boot2])
	PCp_vals2[1:n_boot2,k] = 1-FC_ECDF2[[k]](PCbsample2[k,1:n_boot2])
	PCp_vals3[1:n_boot2,k] = 1-FC_ECDF3[[k]](PCbsample3[k,1:n_boot2])
	PCp_vals4[1:n_boot2,k] = 1-FC_ECDF4[[k]](PCbsample4[k,1:n_boot2])
	PCp_vals=cbind(PCp_vals1,PCp_vals2,PCp_vals3,PCp_vals4)

# compute the minimum p-value based on the original sample
	p_val_1 =  min( 1- FC_ECDF1[[k]](Tn1[k]), p_val_1 )
	p_val_2 =  min( 1- FC_ECDF2[[k]](Tn2[k]), p_val_2 )
	p_val_3 =  min( 1- FC_ECDF3[[k]](Tn3[k]), p_val_3 )
	p_val_4 =  min( 1- FC_ECDF4[[k]](Tn4[k]), p_val_4 )
	p_val =  min(p_val_1,p_val_2,p_val_3,p_val_4)
}  # end of loop of second bootstrap

# compute adjusted p-values - full recentering 
min_p_vec <- apply(FCp_vals,1,min) # minimum p-value vector
edfF_n <- ecdf(min_p_vec)    # empirical distribution function
FCpval = edfF_n(p_val)  

# compute adjusted p-values - partial recentering 
min_p_vec <- apply(PCp_vals,1,min)
edfF_n <- ecdf(min_p_vec)
PCpval = edfF_n(p_val) 

# CHEN AND SZROETER TEST
J=c(n*varTn1b,n*varTn2b,n*varTn3b, n*varTn4b)
J[J==0]=0.000000000001
theta=1/sqrt(J)
#Tuning parameter for the sample size
K=1/delta_n
# Smooth the indicator function
psi=pnorm(K*theta*Tn);
Lambda=(dnorm(K*theta*Tn)*K)/sqrt(n)
# Create a diagonal matrix with the variances of the elements of Tn
delta=diag(theta);
# Compute the P-value
Q1=sqrt(n)*t(psi)%*%(delta%*%Tn)-matrix(1,1,length(Tn))%*%(Lambda);
Q2=sqrt(t(psi)%*%delta%*%t(t(J)%*%delta*psi))
Pc=min(1-pnorm(Q1/max(Q2,0.00000001)),1)

list(ChSz.pval=Pc, fullr.pval=FCpval, partialr.pval=PCpval)

}


# IV test based on Bennett (2009) and Chen and Schroeter (2012)
ivtest<-function(y,d,z,n_boot1=999, n_boot2=999){
n=length(y)
p1_1<-mean(d[z==1])
p1_0<-mean(d[z==0])
p0_1<-1-mean(d[z==1])
p0_0<-1-mean(d[z==0])
z1<-z[d==1]
z0<-z[d==0]
y1<-y[d==1]
y0<-y[d==0]
y11=sort(y1[z1==1])
y01=sort(y1[z1==0])
y10=sort(y0[z0==1])
y00=sort(y0[z0==0])

y11temp<-sort(runif(length(y11)))
y00temp<-sort(runif(length(y00)))

y11min<-(y11[y11temp<=quantile(y11temp,p1_0/p1_1)]) 
y11max<-(y11[y11temp>=quantile(y11temp, 1-p1_0/p1_1)])
y00min<-(y00[y00temp<=quantile(y00temp, p0_1/p0_0)])
y00max<-(y00[y00temp>=quantile(y00temp,1-p0_1/p0_0)])

#compute the normalized differences
normdif.treated<-((mean(y11)>mean(y01))*(mean(y11min)-mean(y01))-(mean(y11)<mean(y01))*(mean(y11max)-mean(y01)))/sd(y)   
normdif.nontreated<-((mean(y00)>mean(y10))*(mean(y00min)-mean(y10))-(mean(y00)<mean(y10))*(mean(y00max)-mean(y10)))/sd(y)

#first bootstrap 
Tn1b = c()
Tn2b = c()
Tn3b = c()
Tn4b = c()

while(length(Tn1b)<n_boot1){
sboot<-sample(1:length(y),length(y),TRUE)
zb<-z[sboot]
db<-d[sboot]
yb<-y[sboot]

obsb=length(yb)
obs1b<-sum(db)
obs0b<-sum(1-db)

p1_1b<-mean(db[zb==1])
p1_0b<-mean(db[zb==0])
p0_1b<-1-mean(db[zb==1])
p0_0b<-1-mean(db[zb==0])

if (p1_0b/p1_1b<=1){
z1b<-zb[db==1]
z0b<-zb[db==0]
y1b<-yb[db==1]
y0b<-yb[db==0]

y11b=sort(y1b[z1b==1])
y01b=sort(y1b[z1b==0])
y10b=sort(y0b[z0b==1])
y00b=sort(y0b[z0b==0])

y11tempb<-sort(runif(length(y11b)))
y00tempb<-sort(runif(length(y00b)))

y11minb<-mean(y11b[y11tempb<=quantile(y11tempb,p1_0b/p1_1b)]) 
y11maxb<-mean(y11b[y11tempb>=quantile(y11tempb, 1-p1_0b/p1_1b)])
y00minb<-mean(y00b[y00tempb<=quantile(y00tempb, p0_1b/p0_0b)])
y00maxb<-mean(y00b[y00tempb>=quantile(y00tempb,1-p0_1b/p0_0b)])

y01b=mean(y01b)
y10b=mean(y10b)

Tn1b = c(Tn1b,(y11minb-y01b)) 
Tn2b = c(Tn2b,(y01b-y11maxb)) 
Tn3b = c(Tn3b,(y00minb-y10b)) 
Tn4b = c(Tn4b,(y10b-y00maxb)) 
}
}

# variances
varTn1b=var(Tn1b)
varTn2b=var(Tn2b)
varTn3b=var(Tn3b)
varTn4b=var(Tn4b)

# create objects for partial and full recentering
FCbsample1 <-matrix(0,nrow=1,ncol=n_boot2)
FCbsample2 <-matrix(0,nrow=1,ncol=n_boot2)
FCbsample3 <-matrix(0,nrow=1,ncol=n_boot2)
FCbsample4 <-matrix(0,nrow=1,ncol=n_boot2)

PCbsample1 <-matrix(0,nrow=1,ncol=n_boot2)
PCbsample2 <-matrix(0,nrow=1,ncol=n_boot2)
PCbsample3 <-matrix(0,nrow=1,ncol=n_boot2)
PCbsample4 <-matrix(0,nrow=1,ncol=n_boot2)

FCp_vals1 <- matrix(0, nrow = n_boot2, ncol = 1)
FCp_vals2 <- matrix(0, nrow = n_boot2, ncol = 1)
FCp_vals3 <- matrix(0, nrow = n_boot2, ncol = 1)
FCp_vals4 <- matrix(0, nrow = n_boot2, ncol = 1)

PCp_vals1 <- matrix(0, nrow = n_boot2, ncol = 1)
PCp_vals2 <- matrix(0, nrow = n_boot2, ncol = 1)
PCp_vals3 <- matrix(0, nrow = n_boot2, ncol = 1)
PCp_vals4 <- matrix(0, nrow = n_boot2, ncol = 1)

#Fix value of sequence delta_n
delta_n = sqrt(2*log(log(n))/n)

# compute test statistics
Tn1 = (mean(y11min)-mean(y01)) 
Tn2 = (mean(y01)-mean(y11max)) 
Tn3 = (mean(y00min)-mean(y10)) 
Tn4 = (mean(y10)-mean(y00max))
Tn=c(Tn1,Tn2, Tn3, Tn4) 

# compute recentering terms
maxVec1=( Tn1 >= (-delta_n*sqrt(varTn1b))  )*Tn1    + ( Tn1 <  (-delta_n*sqrt(varTn1b))  )* (-delta_n*sqrt(varTn1b) )
maxVec2=( Tn2 >= (-delta_n*sqrt(varTn2b))  )*Tn2   +( Tn2 < (-delta_n*sqrt(varTn2b))  )* (-delta_n*sqrt(varTn2b) )
maxVec3=( Tn3 >= (-delta_n*sqrt(varTn3b))  )*Tn3    + ( Tn3 <  (-delta_n*sqrt(varTn3b))  )* (-delta_n*sqrt(varTn3b) )
maxVec4=( Tn4 >= (-delta_n*sqrt(varTn4b))  )*Tn4   +( Tn4 < (-delta_n*sqrt(varTn4b))  )* (-delta_n*sqrt(varTn4b) )

Tn1 = sqrt(n)*Tn1 
Tn2 = sqrt(n)*Tn2
Tn3 = sqrt(n)*Tn3 
Tn4 = sqrt(n)*Tn4

Tn1b=sqrt(n)*Tn1b
Tn2b=sqrt(n)*Tn2b
Tn3b=sqrt(n)*Tn3b
Tn4b=sqrt(n)*Tn4b

# Generate bootstrap samples (using first bootstrap)- full recentering
FCBootstrap1<-matrix(Tn1b-rep(Tn1,n_boot1),1,n_boot1)
FCBootstrap2<-matrix(Tn2b-rep(Tn2,n_boot1),1,n_boot1)
FCBootstrap3<-matrix(Tn3b-rep(Tn3,n_boot1),1,n_boot1)
FCBootstrap4<-matrix(Tn4b-rep(Tn4,n_boot1),1,n_boot1)

# Generate bootstrap samples (using first bootstrap)- partial recentering
PCBootstrap1<-matrix(Tn1b-rep((sqrt(n)* maxVec1),n_boot1),1,n_boot1)
PCBootstrap2<-matrix(Tn2b-rep((sqrt(n)* maxVec2),n_boot1),1,n_boot1)
PCBootstrap3<-matrix(Tn3b-rep((sqrt(n)* maxVec3),n_boot1),1,n_boot1)
PCBootstrap4<-matrix(Tn4b-rep((sqrt(n)* maxVec4),n_boot1),1,n_boot1)

# Generate empirical distribution functions - full recentering
FC_ECDF1 <- apply(FCBootstrap1, 1, ecdf)
FC_ECDF2 <- apply(FCBootstrap2, 1, ecdf)
FC_ECDF3 <- apply(FCBootstrap3, 1, ecdf)
FC_ECDF4 <- apply(FCBootstrap4, 1, ecdf)

# Generate empirical distribution functions - partial recentering
PC_ECDF1 <- apply(PCBootstrap1, 1, ecdf)
PC_ECDF2 <- apply(PCBootstrap2, 1, ecdf)
PC_ECDF3 <- apply(PCBootstrap3, 1, ecdf)
PC_ECDF4 <- apply(PCBootstrap4, 1, ecdf)

# evaluate the test statistics at the empirical distributions to compute the p-values
p_val_1 =  1- FC_ECDF1[[1]](Tn1)
p_val_2 =  1- FC_ECDF2[[1]](Tn2)
p_val_3 =  1- FC_ECDF3[[1]](Tn3)
p_val_4 =  1- FC_ECDF4[[1]](Tn4) 

# compute the minimum p-value based on the original sample
p_val = min(p_val_1,p_val_2,p_val_3,p_val_4)

# second bootstrap - draw observations from first bootstrap
bsample2 <- sample(n_boot1, n_boot2, replace=TRUE)

# draw samples of recentered test statistics out of first bootstrap - full recentering
FCbsample1 = matrix(FCBootstrap1[1,bsample2],1,n_boot2)
FCbsample2 = matrix(FCBootstrap2[1,bsample2],1,n_boot2)
FCbsample3 = matrix(FCBootstrap3[1,bsample2],1,n_boot2)
FCbsample4 = matrix(FCBootstrap4[1,bsample2],1,n_boot2)

# draw samples of recentered test statistics out of first bootstrap - partial recentering
PCbsample1 = matrix(PCBootstrap1[1,bsample2],1,n_boot2)
PCbsample2 = matrix(PCBootstrap2[1,bsample2],1,n_boot2)
PCbsample3 = matrix(PCBootstrap3[1,bsample2],1,n_boot2)
PCbsample4 = matrix(PCBootstrap4[1,bsample2],1,n_boot2)

# empirical distributions of fully recentered p-values
FCp_vals1[1:n_boot2,1] = 1-FC_ECDF1[[1]](FCbsample1[1,1:n_boot2])	
FCp_vals2[1:n_boot2,1] = 1-FC_ECDF2[[1]](FCbsample2[1,1:n_boot2])
FCp_vals3[1:n_boot2,1] = 1-FC_ECDF3[[1]](FCbsample3[1,1:n_boot2])
FCp_vals4[1:n_boot2,1] = 1-FC_ECDF4[[1]](FCbsample4[1,1:n_boot2])
FCp_vals=cbind(FCp_vals1,FCp_vals2,FCp_vals3,FCp_vals4)

# empirical distributions of partially recentered p-values
PCp_vals1[1:n_boot2,1] = 1-FC_ECDF1[[1]](PCbsample1[1,1:n_boot2])
PCp_vals2[1:n_boot2,1] = 1-FC_ECDF2[[1]](PCbsample2[1,1:n_boot2])
PCp_vals3[1:n_boot2,1] = 1-FC_ECDF3[[1]](PCbsample3[1,1:n_boot2])
PCp_vals4[1:n_boot2,1] = 1-FC_ECDF4[[1]](PCbsample4[1,1:n_boot2])
PCp_vals=cbind(PCp_vals1,PCp_vals2,PCp_vals3,PCp_vals4)

# compute adjusted p-values - full recentering 
min_p_vec <- apply(FCp_vals,1,min) # minimum p-value vector
edfF_n <- ecdf(min_p_vec)    # empirical distribution function
FCpval = edfF_n(p_val)  

# compute adjusted p-values - partial recentering 
min_p_vec <- apply(PCp_vals,1,min)
edfF_n <- ecdf(min_p_vec)
PCpval = edfF_n(p_val) 

# CHEN AND SZROETER TEST
J=c(n*varTn1b,n*varTn2b,n*varTn3b, n*varTn4b)
J[J==0]=0.000000000001
theta=1/sqrt(J)
#Tuning parameter for the sample size
K=1/delta_n
# Smooth the indicator function
psi=pnorm(K*theta*Tn);
Lambda=(dnorm(K*theta*Tn)*K)/sqrt(n)
# Create a diagonal matrix with the variances of the elements of Tn
delta=diag(theta);
# Compute the P-value
Q1=sqrt(n)*t(psi)%*%(delta%*%Tn)-matrix(1,1,length(Tn))%*%(Lambda);
Q2=sqrt(t(psi)%*%delta%*%t(t(J)%*%delta*psi))
Pc=min(1-pnorm(Q1/max(Q2,0.00000001)),1)

list(ChSz.pval=Pc, fullr.pval=FCpval, partialr.pval=PCpval, normdif.treated=normdif.treated, normdif.nontreated=normdif.nontreated)
} 


# IV test based on Bonferroni adjustment
ivtest.bonferroni<-function(y,d,z,m, boot=499){
obs=length(y)
obs1<-sum(d)
obs0<-sum(1-d)
p1_1<-mean(d[z==1])
p1_0<-mean(d[z==0])
p0_1<-1-mean(d[z==1])
p0_0<-1-mean(d[z==0])
z1<-z[d==1]
z0<-z[d==0]
y1<-y[d==1]
y0<-y[d==0]
y11=sort(y1[z1==1])
y01=sort(y1[z1==0])
y10=sort(y0[z0==1])
y00=sort(y0[z0==0])

y11temp<-sort(runif(length(y11)))
y00temp<-sort(runif(length(y00)))

y11min<-(y11[y11temp<=quantile(y11temp,p1_0/p1_1)]) 
y11max<-(y11[y11temp>=quantile(y11temp, 1-p1_0/p1_1)])
y00min<-(y00[y00temp<=quantile(y00temp, p0_1/p0_0)])
y00max<-(y00[y00temp>=quantile(y00temp,1-p0_1/p0_0)])

stat1<- (mean(y11min)-mean(y01))
stat2<- (mean(y01)-mean(y11max))
stat3<- (mean(y00min)-mean(y10))
stat4<- (mean(y10)-mean(y00max))

stat1b<-c()
stat2b<-c()
stat3b<-c()
stat4b<-c()

while(length(stat1b)<boot){
sboot<-sample(1:obs,m,TRUE)
zb<-z[sboot]
db<-d[sboot]
yb<-y[sboot]

obsb=length(yb)
obs1b<-sum(db)
obs0b<-sum(1-db)

p1_1b<-mean(db[zb==1])
p1_0b<-mean(db[zb==0])
p0_1b<-1-mean(db[zb==1])
p0_0b<-1-mean(db[zb==0])

if ((p1_0b/p1_1b<=1) & (p0_1b/p0_0b<=1)){
z1b<-zb[db==1]
z0b<-zb[db==0]
y1b<-yb[db==1]
y0b<-yb[db==0]
y11b=sort(y1b[z1b==1])
y01b=sort(y1b[z1b==0])
y10b=sort(y0b[z0b==1])
y00b=sort(y0b[z0b==0])

y11tempb<-sort(runif(length(y11b)))
y00tempb<-sort(runif(length(y00b)))

y11minb<-mean(y11b[y11tempb<=quantile(y11tempb,p1_0b/p1_1b)]) 
y11maxb<-mean(y11b[y11tempb>=quantile(y11tempb, 1-p1_0b/p1_1b)])
y00minb<-mean(y00b[y00tempb<=quantile(y00tempb, p0_1b/p0_0b)])
y00maxb<-mean(y00b[y00tempb>=quantile(y00tempb,1-p0_1b/p0_0b)])

stat1b<-c(stat1b, (mean(y11minb)-mean(y01b))-stat1)
stat2b<-c(stat2b, (mean(y01b)-mean(y11maxb))-stat2)
stat3b<- c(stat3b, (mean(y00minb)-mean(y10b))-stat3)
stat4b<- c(stat4b, (mean(y10b)-mean(y00maxb))-stat4)
}
}
temp<-is.na(stat1b)
stat1b<-stat1b[temp==0]
temp<-is.na(stat2b)
stat2b<-stat2b[temp==0]
temp<-is.na(stat3b)
stat3b<-stat3b[temp==0]
temp<-is.na(stat4b)
stat4b<-stat4b[temp==0]

sd1<-sd(stat1b)
sd2<-sd(stat2b)
sd3<-sd(stat3b)
sd4<-sd(stat4b)

normdif.treated=max(stat1/sd(y),stat2/sd(y))
normdif.nontreated=max(stat3/sd(y),stat4/sd(y))

stat1b<-stat1b/sd1
stat2b<-stat2b/sd2
stat3b<-stat3b/sd3
stat4b<-stat4b/sd4

stat1<-stat1/sd1
stat2<-stat2/sd2
stat3<-stat3/sd3
stat4<-stat4/sd4

pval1=(stat1>0)*(1-mean(stat1>stat1b))+(stat1<=0)
pval2=(stat2>0)*(1-mean(stat2>stat2b))+(stat2<=0)
pval3=(stat3>0)*(1-mean(stat3>stat3b))+(stat3<=0)
pval4=(stat4>0)*(1-mean(stat4>stat4b))+(stat4<=0)
pval=min(pval1,pval2,pval3,pval4)
pval=4*pval
list(pval.adj=pval, normdif.treated=normdif.treated, normdif.nontreated=normdif.nontreated)
}


# LATE estimator
late<-function(y,d,z){
late<-(mean(y[z==1])-mean(y[z==0]))/(mean(d[z==1])-mean(d[z==0]))
list (late=late)
}

latenomon<-function(y,d,z,discrete=FALSE, bwp1=NULL, bwp0=NULL, bwq1=NULL, bwq0=NULL, multiplier=1,nn=length(y)){
ymax<-max(y)
ymin<-min(y)
p1_1<-mean(d[z==1])
p1_0<-mean(d[z==0])
d1=d[z==1]
d0=d[z==0]
y1=y[z==1]
y0=y[z==0]
p1<-c()
p0<-c()
q1<-c()
q0<-c()
maxp1q1<-c()
minp1q1<-c()
maxp0q0<-c()
minp0q0<-c()

if (discrete==1){
masspoints<-seq(from=ymin,to=ymax,by=1)
for (i in 1:length(masspoints)){
p1= c(p1,mean( (y1==masspoints[i])& d1==1))
p0= c(p0,mean( (y1==masspoints[i])& d1==0))
q1= c(q1,mean( (y0==masspoints[i])& d0==1))
q0= c(q0,mean( (y0==masspoints[i])& d0==0))
maxp1q1<-c(maxp1q1,max(p1[i],q1[i]))
minp1q1<-c(minp1q1,min(p1[i],q1[i]))
maxp0q0<-c(maxp0q0,max(p0[i],q0[i]))
minp0q0<-c(minp0q0,min(p0[i],q0[i]))
}
lambda1=sum(minp1q1)
lambda0=sum(minp0q0)
delta0=1-sum(maxp0q0)
delta1=1-sum(maxp1q1)
late.cd<-(sum(masspoints*(maxp1q1-minp1q1))-sum(masspoints*(maxp0q0-minp0q0)))/(p1_1+p1_0-2*lambda1)
late.c<-(sum(masspoints*(p1-minp1q1))-sum(masspoints*(q0-minp0q0)))/(p1_1-lambda1)
late.d<-(sum(masspoints*(q1-minp1q1))-sum(masspoints*(p0-minp0q0)))/(p1_0-lambda1)
late.mon<-(mean(y[z==1])-mean(y[z==0]))/(mean(d[z==1])-mean(d[z==0]))
}
else{
temp<-is.null(bwp1)
if (temp==1){
bwp1 <- npudensbw(formula=~y1[d1==1],bwmethod="normal-reference")
}
#p1 <-npudens(bws=bwp1, edat=seq(ymin,ymax,length.out = n))$dens
p1<-density(x=y1[d1==1], bw = bwp1$bw, n = nn, from=ymin, to=ymax, adjust = multiplier)$y
p1 <- p1/sum(p1)*p1_1

#p1<-density(x=y1[d1==1], bw = bwp1$bw, n = n, from=ymin, to=ymax)$y


temp<-is.null(bwp0)
if (temp==1){
bwp0 <- npudensbw(formula=~y1[d1==0],bwmethod="normal-reference")
}
#p0 <-npudens(bws=bwp0, edat=seq(ymin,ymax,length.out = n))$dens
p0<-density(x=y1[d1==0], bw = bwp0$bw, n = nn, from=ymin, to=ymax, adjust = multiplier)$y
#p0 <-npudens(bws=bwp0,edat=y)$dens
p0 <- p0/sum(p0)*(1-p1_1)

temp<-is.null(bwq1)
if (temp==1){
bwq1 <- npudensbw(formula=~y0[d0==1], bwmethod="normal-reference")
}
#q1 <-npudens(bws=bwq1, edat=seq(ymin,ymax,length.out = n))$dens
q1<-density(x=y0[d0==1], bw = bwq1$bw, n = nn, from=ymin, to=ymax, adjust = multiplier)$y
#q1 <-npudens(bws=bwq1, edat=y)$dens
q1 <- q1/sum(q1)*p1_0

temp<-is.null(bwq0)
if (temp==1){
bwq0 <- npudensbw(formula=~y0[d0==0],bwmethod="normal-reference")
}
#q0 <-npudens(bws=bwq0, edat=seq(ymin,ymax,length.out = n))$dens
q0<-density(x=y0[d0==0], bw = bwq0$bw, n = nn, from=ymin, to=ymax, adjust = multiplier)$y
#q0 <-npudens(bws=bwq0,edat=y)$dens
q0 <- q0/sum(q0)*(1-p1_0)

for (i in 1:nn){
maxp1q1<-c(maxp1q1,max(p1[i],q1[i]))
minp1q1<-c(minp1q1,min(p1[i],q1[i]))
maxp0q0<-c(maxp0q0,max(p0[i],q0[i]))
minp0q0<-c(minp0q0,min(p0[i],q0[i]))
}
lambda1=sum(minp1q1)
lambda0=sum(minp0q0)
delta0=1-sum(maxp0q0)
delta1=1-sum(maxp1q1)
late.mon<-(mean(y[z==1])-mean(y[z==0]))/(mean(d[z==1])-mean(d[z==0]))
y=seq(ymin,ymax,length.out=nn)
late.cd<-(sum(y*(maxp1q1-minp1q1))-sum(y*(maxp0q0-minp0q0)))/(p1_1+p1_0-2*lambda1)
late.c<-(sum(y*(p1-minp1q1))-sum(y*(q0-minp0q0)))/(p1_1-lambda1)
late.d<-(sum(y*(q1-minp1q1))-sum(y*(p0-minp0q0)))/(p1_0-lambda1)
}
list(lambda1=lambda1, lambda0=lambda0, delta0=delta0, late.cd=late.cd, late.c=late.c, late.d=late.d, late.mon=late.mon, c1=p1_1-min(lambda1,delta0), c0=1-p1_0-min(lambda0,delta1), d1=sum(q1-minp1q1), d0=sum(p0-minp0q0), bwp1=bwp1, bwp0=bwp0, bwq1=bwq1, bwq0=bwq0)  
}

bootstrap.late<-function(y,d,z,boot=1999,discrete=FALSE, bwp1, bwp0, bwq1, bwq0, multiplier=1, nn=length(y)){
obs<-length(y)
mc=c()
temp=c()
while(length(temp)<boot){
sboot<-sample(1:obs,obs,TRUE)
zb<-z[sboot]
db<-d[sboot]
yb=y[sboot]
mc<-rbind(mc,c(latenomon(y=yb,d=db,z=zb,discrete=discrete, bwp1=bwp1, bwp0=bwp0, bwq1=bwq1, bwq0=bwq0, multiplier=multiplier,nn=nn)))
temp<-c(temp,1)
}
list(mc=mc, se.lambda1=sd(as.numeric(mc[,1])), se.delta0=sd(as.numeric(mc[,2])), se.late.cd=sd(as.numeric(mc[,3])), se.late.c=sd(as.numeric(mc[,4])), se.late.d=sd(as.numeric(mc[,5])), se.late.mon=sd(as.numeric(mc[,6])), se.c1=sd(as.numeric(mc[,7])), se.c0=sd(as.numeric(mc[,8])), se.d1=sd(as.numeric(mc[,9])), se.d0=sd(as.numeric(mc[,10])), se.diff.lambda1.delta0=sd(as.numeric(mc[,1])-as.numeric(mc[,2])) )  
}

bootstrap.late.sub<-function(y,d,z,boot=1999,discrete=FALSE, bwp1=NULL, bwp0=NULL, bwq1=NULL, bwq0=NULL, multiplier=1,m, nn=length(y)){
mc=c()
temp=c()
while(length(temp)<boot){
sboot<-sample(1:length(y),m,FALSE)
zb<-z[sboot]
db<-d[sboot]
yb=y[sboot]
mc<-rbind(mc,c(latenomon(yb,db,zb,discrete=discrete, bwp1=bwp1, bwp0=bwp0, bwq1=bwq1, bwq0=bwq0, multiplier=multiplier, nn=nn)))
temp<-c(temp,1)
}

list(mc=mc, se.lambda1=sd(as.numeric(mc[,1])), se.one_delta0=sd(as.numeric(mc[,2])), se.late.cd=sd(as.numeric(mc[,3])), se.late.c=sd(as.numeric(mc[,4])), se.late.d=sd(as.numeric(mc[,5])), se.late.mon=sd(as.numeric(mc[,6])), se.c1=sd(as.numeric(mc[,7])), se.c0=sd(as.numeric(mc[,8])), se.d1=sd(as.numeric(mc[,9])), se.d0=sd(as.numeric(mc[,10])), se.diff.lambda1.one_delta0=sd(as.numeric(mc[,7])-as.numeric(mc[,8])) )  
}


###############
# Application #
###############


data=read.dta("C:\\Card IV data.dta")
attach(data)

z=nearc4
d=degree
y=lwage

#testing inequality constraints
ivtest.bonferroni(y,d,z, m=length(y), boot=1999)
ivtestprob(y,d,z, n_boot1=1999, n_boot2=1999, subsets=2)
ivtest(y,d,z, n_boot1=1999,  n_boot2=1999)

# testing equality constraints
z1<-z[d==1]
z0<-z[d==0]
y1<-y[d==1]
y0<-y[d==0]
y11=y1[z1==1]
y01=y1[z1==0]
y10=y0[z0==1]
y00=y0[z0==0]
t.test(x=y11, y=y01)
t.test(x=y00, y=y10)

# testing in subsamples
indicator=(fatheduc<12 &  black==0 &  south66==0 &  smsa66==1 &  misfath==0)
z=nearc4[indicator==1]
d=degree[indicator==1]
y=lwage[indicator==1]
ivtest(y,d,z, n_boot1=1999,  n_boot2=1999)




##############
# Simulation #
##############

mctestpartial<-c( )
mctestfull<-c( )
mcChenSzroeter<-c( )
mcdiff1<-c( )
mcdiff2<-c( )
mclate<-c()
mcnaive<-c()
mccomplier<-c()

n<-1000
sigma <- matrix(c(1,0.5,0.5,1),2,2)
repetitions<-1000
i=1
while(length(mctestpartial)<repetitions){
set.seed(i)
e<-(rmvnorm(n,rep(0,2),sigma))
z<-rnorm(n)>0
d<-(0.5*z+e[,1])>0
y<-d+z+e[,2]

p1_1<-mean(d[z==1])
p1_0<-mean(d[z==0])
p0_1<-1-mean(d[z==1])
p0_0<-1-mean(d[z==0])

z1<-z[d==1]
z0<-z[d==0]
y1<-y[d==1]
y0<-y[d==0]
y11=y1[z1==1]
y01=y1[z1==0]
y10=y0[z0==1]
y00=y0[z0==0]
mclate=c(mclate, late(y,d,z)$late)
mcnaive<-c(mcnaive, mean(y[d==1])-mean(y[d==0]))
mccomplier<-c(mccomplier, mean(d[z==1])-mean(d[z==0]))
temp<-ivtest(y,d,z, n_boot1=499, n_boot2=499)
mctestpartial=c(mctestpartial, temp$partialr.pval)
mctestfull=c(mctestfull, temp$fullr.pval)
mcChenSzroeter<-c(mcChenSzroeter,  temp$ChSz.pval)
mcdiff1<-c(mcdiff1, temp$normdif.treated) 
mcdiff2<-c(mcdiff2, temp$normdif.nontreated)
i=i+1
}

mean(mclate)
mean(mcnaive)
mean(mcdiff1)
mean(mcdiff2)
mean(mctestpartial<0.05)
mean(mctestfull<0.05)
mean(mcChenSzroeter<0.05)